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LONG-TERM  GOALS 

The  long-term  goal  of  this  project  is  to  design  and  evaluate  the  components  that  will  comprise  a  next  genera¬ 
tion  mesoscale  atmospheric  model  within  the  Coupled  Ocean/Atmosphere  Mesoscale  Prediction  System 
(COAMPS®1).  It  is  anticipated  that  in  order  to  meet  future  Navy  requirements,  next  generation  approaches  to 
numerical  techniques  and  physical  parameterizations  will  be  needed. 

OBJECTIVES 

The  objectives  of  this  project  involve  the  development,  testing,  and  validation  of:  i)  new  numerical 
techniques  such  as  advection  schemes  and  time  differencing  methods,  and  ii)  new  methods  for  repre¬ 
senting  cloud-scale  physical  processes.  Both  of  these  objectives  are  tailored  to  address  high-resolution 
applications  for  horizontal  grid  increments  at  1  km  or  less. 

APPROACH 

Our  approach  is  to  follow  a  methodical  plan  in  the  development  and  testing  of  a  nonhydro static  micro¬ 
scale  modeling  system  that  will  leverage  the  existing  COAMPS  and  new  model  prototypes.  Our  work 
on  numerical  methods  will  involve  investigation  of  spatial  and  temporal  discretization  algorithms  that 
are  superior  to  the  current  generation  leap-frog,  second-order  accurate  numerical  techniques  presently 
employed  in  COAMPS  and  many  other  models;  these  new  discretization  methods  will  be  developed 
and  implemented.  Our  work  on  the  physics  for  the  next-generation  COAMPS  will  feature  the  devel¬ 
opment  of  physical  parameterizations  specifically  designed  to  represent  cloud-scale  processes  operat¬ 
ing  on  fine  scales.  A  parameterization  is  proposed  that  properly  represents  the  coupled  nature  between 
the  turbulence  and  microphysics  in  droplet  activation,  evaporation,  and  auto-conversion  processes  for 
mesoscale  and  microscale  models.  Validation  and  evaluation  of  the  modeling  system  will  be  per¬ 
formed  using  datasets  of  opportunity,  particularly  in  regions  of  Navy  significance. 

WORK  COMPLETED 

1.  Development  of  a  new  parameterization  of  turbulent  cloud  droplet  number  density  flux. 

Prediction  of  cloud  microphysics  is  a  crucial  part  of  next  generation  COAMPS;  it  depends  on  the  real¬ 
istic  parameterization  of  turbulent  mixing  of  cloud  droplets.  Traditionally,  the  turbulent  mixing  of  the 
droplet  number  concentration  is  calculated  based  on  the  so-called  down-gradient  method  in  which  the 
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droplets  are  diffused  along  the  direction  of  decreasing  concentration.  This  approach,  however,  is  se¬ 
verely  flawed,  because  it  does  not  include  key  physical  processes  that  are  critical  to  the  droplet  spec¬ 
trum  evolution,  i.e.  the  droplet  activation  in  turbulent  updrafts  and  the  evaporation  in  the  downdrafts. 
A  new  mixing  scheme  is  developed  to  include  these  processes;  it  can  be  written  as 


dw'N'  „ 

- =  Transport  +  dissipation  + 


where  w'N'  is  turbulent  flux  of  the  droplet  number  density,  first  two  tenns  on  the  right-hand- side 
represents  normal  processes  related  turbulence  transport  and  dissipation,  and  the  last  is  the  activation 
contribution  to  the  flux.  The  activation  term  is  formulated  in  tenns  of  the  probability  density  function 
of  vertical  velocity  and  aerosol  activation  spectrum.  This  tenn  makes  significant  contribution  at  the 
cloud  base  where  the  activation  occurs  in  the  updrafts.  Therefore,  the  turbulent  flux  is  closely  coupled 
to  the  turbulence  structure.  We  included  this  parameterization  in  a  single  column  high-resolution  tur¬ 
bulence  closure  model  and  tested  the  perfonnance  against  DYCOMS  (Dynamics  and  Chemistry  of 
Marine  Stratocumulus  Cloud  Experiment  II)  observations. 

2.  New  lower  boundary  condition  for  the  COAMPS  dynamical  core 

Over  the  past  several  years,  a  negative  geopotential  height  bias  has  been  apparent  in  the  COAMPS 
model  during  wintertime  months.  This  height  bias  was  particularly  apparent  in  mid-latitude  regions 
such  as  the  FNMOC  operational  Europe  area.  Height  biases  of  -50  m  at  500  hPa  were  common  the 
wintertime  months  and  much  larger  than  the  NOGAPS  height  biases.  Over  the  past  year,  a  series  of 
real  data  and  simplified  tests  were  run  using  CO  AMPS,  the  results  of  which  motivated  the  implemen¬ 
tation  of  a  new  lower  boundary  condition  for  the  COAMPS  nonhydrostatic  model. 

3.  Spectral  element  and  discontinuous  Galerkin  2D  prototypes. 

In  FY07  we  continue  examination  of  three  different  nonhydrostatic  equation  sets  and  explore  solving 
these  equation  sets  using  spectral  element  and  discontinuous  Galerkin  techniques.  We  have  continued 
collaboration  with  Frank  Giraldo  at  NPS  with  this  work.  Spectral  element  (SE)  (see  Giraldo  2005)  and 
discontinuous  Galerkin  (DG)  (see  Giraldo  2006)  methods  are  a  new  class  of  spatial  discretization 
methods  that  are  used  to  approximate  the  derivatives  of  the  governing  equations.  The  advantage  of  SE 
and  DG  methods  is  that  they  offer  high-order  accuracy  and  this  accuracy  can  be  achieved  on  any  un¬ 
structured  grid  -  this  is  not  true  for  finite  difference  methods  where  the  differencing  stencils  require  a 
certain  level  of  structure  (such  as  orthogonality).  Based  on  SE  and  DG  methods,  we  examined  5  differ¬ 
ent  prototypes  for  the  nonhydrostatic  Euler  equations.  We  tested  various  aspects  of  these  models  using 
idealized  cases  for  warm  bubbles,  gravity  currents,  and  gravity  waves. 

4.  Development  of  new  cloud  microphysical  equations  with  variable  collection  efficiencies 

For  many  applications  of  interest,  it  is  necessary  to  consider  mixed-phase  cloud  microphysical  proc¬ 
esses.  In  standard  bulk  cloud  microphysical  schemes,  collection  rates  are  computed  between  the  ice 
and  liquid  condensate  categories  with  the  result  often  being  the  creation  a  third  condensate  category 
(such  as  rain  and  snow  collisions  leading  to  the  formation  of  graupel).  Such  rate  calculations  have 
been  shown  to  lead  to  the  over-depletion  of  the  collected  scalar  within  a  given  time  step  (referred  to  as 
the  over-depletion  problem).  Modelers  typically  resort  to  placing  add  hoc  limits  to  the  rate  equations  to 
preserve  a  positive-definite  solution.  A  recent  study  by  Gaudet  and  Schmidt  (2007)  analyzed  the  na¬ 
ture  of  these  errors  and  derived  a  new  fonnulization  of  the  rate  equations,  termed  “Numerical  Bound- 
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ing”,  which  is  based  on  accurate  numerical  integration  of  the  collection  terms.  The  new  equations  ac¬ 
curately  account  for  errors  in  the  standard  time  discretization  of  the  collection  equations  allowing 
much  larger  time  steps  to  be  taken,  a  fact  which  may  prove  useful  for  envisioned  semi-Lagrangian  nu¬ 
merical  frameworks.  We  have  now  extended  these  results  to  include  variable  collection  efficiencies 
that  account  for  known  drop/crystal  interactions.  The  work  has  involved  developing  a  parameterization 
of  the  collection  efficiency  calculations  of  Wang  and  Ji  (2000)  for  a  variety  of  crystal  and  droplet  di¬ 
ameters. 

5.  Weighted  Essentially  Non-Oscillatory  (WENO)  methods. 

Atmospheric  models  require  numerical  methods  that  can  accurately  represent  the  transport  of  tracers 
with  steep  gradients,  such  as  those  that  occur  at  cloud  boundaries  or  the  edges  of  chemical  plumes.  In 
atmospheric  sciences,  the  most  widely  used  numerical  techniques  for  this  type  of  problem  are  flux- 
corrected  transport  or  closely  related  flux-limiter  methods.  The  limiters  are  typically  designed  to  pre¬ 
vent  the  development  of  new  extrema  in  the  concentration  field.  This  will  preserve  the  non-negativity 
of  initially  non-negative  fields,  which  is  essential  for  the  correct  simulation  of  cloud  microphysics  or 
chemical  reactions.  One  serious  systematic  weakness  of  flux  limiter  methods  is  that  they  also  tend  to 
damp  the  amplitude  of  extrema  in  smooth  regions  of  the  flow,  such  as  the  trough  of  a  well-resolved 
sine  wave.  To  avoid  this  problem,  we  have  been  investigating  the  application  of  WENO  (Weighted 
Essentially  Non-Oscillatory)  methods  to  tracer  transport  in  atmospheric  models.  WENO  methods  are 
widely  used  in  many  disciplines,  but  scarcely  been  tested  in  atmospheric  applications.  WENO  meth¬ 
ods  preserve  steep  gradients  while  simultaneously  avoiding  the  dissipation  of  smooth  extrema  by  esti¬ 
mating  the  value  of  the  solution  in  a  way  that  heavily  weights  the  smoothest  possible  cubic  polynomial 
fit  to  the  local  function  values.  Where  the  solution  is  well  resolved,  all  possible  cubic  interpolants  are 
weighted  almost  equally.  Near  a  steep  gradient,  those  interpolants  are  almost  completely  ignored. 

RESULTS 

1.  Development  of  a  new  parameterization  of  turbulent  cloud  droplet  number  density  flux. 

Fig.  1  shows  the  comparison  of  the  results  between  the  old  and  new  parameterizations.  The  traditional 
“down-gradient”  scheme  predicts  lower  cloud  droplet  number  density  than  the  new  parameterization; 
the  result  predicted  by  the  new  scheme  is  closer  to  the  observation  (circles  with  horizontal  bars).  This 
difference  comes  from  the  different  treatment  in  the  turbulent  droplet  number  flux.  The  old  “down- 
gradient”  approach  diffuses  significant  amount  of  droplets  down  to  the  cloud  base  where  they  are 
evaporated  and  exit  clouds;  the  new  scheme  produces  “up-gradient”  flux  profile  because  of  the  activa¬ 
tion  in  the  updrafts,  leading  to  upward  transport  of  the  droplets  within  cloud  layer.  It  is  seen  that  the 
new  flux  scheme  indeed  improves  the  cloud  microphysical  prediction.  The  parameterized  turbulent 
flux  is  also  compared  with  that  computed  from  NRL  COAMPS-LES  model.  Since  a  LES  model  can 
explicitly  resolve  the  stochastic  motion  of  turbulent  eddies,  it  can  correctly  simulate  the  droplet  num¬ 
ber  flux.  Thus  COAMPS-LES  provides  a  basic  evaluation  of  the  parameterization  result.  It  is  seen  in 
Fig.  2  that  the  parameterized  turbulent  flux  profile  is  very  close  to  that  from  COAMPS-LES.  Both 
COAMPS-LES  and  the  new  parameterization  predict  the  “up-gradient”  feature  just  above  the  cloud 
base,  which  further  confirms  our  view  that  the  new  method  is  more  realistic  than  the  old  down-gradient 
approach. 

2.  New  lower  boundary  condition  for  the  COAMPS  dynamical  core 

During  FY07,  a  series  of  real  data  and  simplified  tests  were  run  using  COAMPS  in  order  to  determine 
the  source  of  a  geopotential  height  bias.  After  a  series  of  intensive  sensitivity  tests,  it  was  determined 
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that  the  formulation  of  the  lower  boundary  condition  for  the  vertical  velocity  equation  was  the  cause  of 
the  geopotential  height  bias.  Three  different  formulations  were  tested  for  the  lower  boundary  condi¬ 
tion.  A  second  order  finite  difference  technique  was  found  to  produce  the  best  results  with  the  lowest 
height  bias  statistics.  For  example,  the  mean  sea  level  pressure  bias,  shown  in  Fig.  3,  for  the  period  1- 
14  February  2007  indicates  that  the  new  lower  boundary  condition  reduces  the  sea  level  pressure  bias 
by  approximately  50%  for  the  Europe  operational  area.  The  geopotential  height  bias  and  RMS  errors 
in  the  troposphere  were  reduced  considerably  by  50%  or  more. 

3.  Spectral  element  2D  prototypes. 

In  order  to  determine  which  equation  set  along  with  which  numerical  method  is  best  suited  for  building 
a  next-generation  mesoscale  model,  we  ran  simulations  for  a  number  of  standard  test  cases.  Here  we 
illustrate  one  of  these  test  cases  for  linear  hydrostatic  waves  using  the  spectral  element  dynamical  core 
with  equation  set.  The  setup  of  the  computational  grid  was  modified  to  allow  for  an  arbitrary  stretch¬ 
ing  of  grid  elements  in  both  vertical  and  horizontal  direction.  Consecutively,  the  spatial  resolution  (ei¬ 
ther  vertical  or  horizontal)  is  refined  in  elements  with  smaller  extent.  In  the  example  shown  (Fig.  4), 
the  stretching  was  applied  in  both  horizontal  and  vertical  direction  and  the  grid  elements  are  shallower 
closer  to  the  bottom  boundary,  and  narrower  closer  to  the  obstacle  at  the  center  of  the  domain.  The 
finer  spatial  resolution  in  areas  of  interest  (e.g.  boundary  layer,  lee  slope  of  the  mountain)  also  puts  an 
additional  restriction  on  the  time  step  of  the  model,  on  top  of  the  temporal  restriction  due  to  irregularly 
spaced  nodal  points  within  each  element.  In  this  example,  the  polynomial  order  was  8  in  both  direc¬ 
tions.  The  simulation  depicts  a  linear  hydrostatic  mountain  wave  case  in  an  isentropic  atmosphere, 
with  a  small  obstacle  (height  =  1  m,  half- width  =  10000  m)  at  the  center  of  the  domain.  The  vertical 
velocity,  shown  in  Fig.  5,  agrees  with  analytic  solutions,  even  with  the  coordinate  stretching. 

4.  Development  of  new  cloud  microphysical  equations  with  variable  collection  efficiencies 
Figure  6  shows  the  collision  efficiency  as  computed  by  Wang  and  Ji  (2000)  together  with  the  param¬ 
eterized  efficiency  for  a  crystal  radius  of  800  microns.  The  plots  show  the  rapid  increase  of  the 
drop/ice  collision  efficiency  (ESR)  to  a  value  near  1.0  as  the  droplet  size  increases  beyond  about  30 
microns.  The  ESR  values  remain  steady  near  the  maximum  until  the  drop  fall  speed  approaches  that  of 
the  given  crystal  at  which  point  the  ESR  values  rapidly  fall  to  zero.  An  inspection  of  Fig.  6a  reveals 
that  ESR  remains  fixed  at  a  value  of  zero  for  all  droplets  greater  than  the  critical  diameter.  A  study  of 
drop-drop  interactions  by  Pinsky  et.  al.  2001  suggests  that  the  collision  efficiencies  may  rapidly  in¬ 
crease  again  beyond  this  critical  threshold  to  a  value  near  one  as  the  larger  faster  moving  drops  collide 
with  the  smaller  drops.  We  thus  seek  to  test  the  impact  of  a  similar  effect  in  our  drop/crystal  collision 
efficiency  parameterization  caused  by  faster  moving  droplets  impacting  slower  moving  crystals.  The 
effects  of  applying  a  variable  ESR  to  the  case  of  rain  collecting  snow  are  shown  in  Fig.  7.  The  Figure 
shows  that  by  far  the  greatest  over-depletion  arises  from  the  standard  analytical  solution  (ESR  =  1.0) 
as  it  appears  in  our  operational  COAMPS  code  (solid  line  labeled  A).  By  invoking  a  numerical  solu¬ 
tion  of  the  rain/snow  collection  while  retaining  a  fixed  ESR  =  1 .0,  we  find  that  the  over-depletion  is 
effectively  eliminated  for  a  wide  range  of  rain  mixing  values  (curve  B).  This  reduction  results  in  large 
part  from  the  removal  of  the  so-called  Wisner  approximation  during  the  numerical  integration  (in  the 
analytical  solution,  the  velocity  differences  between  each  species  are  assumed  constant  and  removed 
from  the  integration).  Including  the  numerical  bounding  technique  (yet  keeping  a  fixed  ESR  =  1.0) 
completely  eliminates  the  over-depletion  as  demonstrated  by  Gaudet  and  Schmidt  (2007)  [curve  C]. 
Finally,  adding  a  variable  ESR  as  depicted  in  Fig.  7  to  the  numerically  bounded  solution  results  in  a 
further  reduction  of  this  particular  rate  calculation  to  a  value  quite  near  zero  (curve  D).  Comparing  D 
with  A  we  find  well  over  an  order  of  magnitude  reduction  in  this  calculation  (a  known  graupel  source 
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term).  By  allowing  ESR  to  return  to  a  value  of  one  beyond  the  critical  drop  radius,  however,  leads  to 
nearly  an  identical  solution  to  that  of  curve  C.  This  suggests  that  a  variable  ESR  will  have  its  greatest 
impact  on  snow  or  pristine  ice  collisions  with  cloud  droplet  due  to  the  slower  droplet  fallspeeds. 

5.  Weighted  Essentially  Non-Oscillatory  (WENO)  methods 

We  have  developing  a  new  method,  in  which  WENO  like  criteria  for  the  presence  of  poorly  resolved 
steep  gradients  are  evaluated  in  a  highly  efficient  manner  and  used  to  detennine  whether  spurious  new 
maxima  and  minima  are  likely  to  be  created  in  the  vicinity  of  a  steep  gradient.  This  new  approach, 
which  we  call  "selective  limiting,"  can  be  applied  to  both  flux-corrected  transport  (FCT)  methodolo¬ 
gies  or  to  piecewise  polynomial  approximations  for  the  fluxes  in  finite-volume  methods,  particularly 
the  piecewise  parabolic  (PPM)  or  cubic  methods  (PCM).  Over  a  wide  range  of  tests,  these  schemes 
were  found  to  give  more  accurate  solutions  than  WENO  methods  for  considerably  less  computational 
cost. 

Figure  8  shows  the  basic  behavior  of  the  exact  solution.  An  initially  round  concentration  of  chemical 
tracer  is  deformed  into  a  crescent  in  a  sheared  circular  flow.  At  t=0,  the  flow  in  the  outer  part  of  the 
domain  is  faster,  and  the  flow  near  the  center  is  slower,  than  that  associated  with  solid  body  rotation. 
However  the  strength  of  the  sheared  component  of  the  flow  periodically  reverses  to  that  at  every  mul¬ 
tiple  of  T/2  time  units,  the  total  circular  arc  traveled  by  each  fluid  parcel  is  identical  to  that  for  one 
revolution  of  solid  body  rotation,  and  the  correct  tracer  field  distribution  is  exactly  the  same  as  the  ini¬ 
tial  condition.  Figure  9  compares  numerical  results  generated  by:  (a)  the  unlimited  PCM  method,  (b) 
the  selectively  limited  PCM  method  with  positivity  preservation,  (d)  the  standard  WENO  method,  and 
(e)  the  selectively  limited  PCM  method  using  a  semi-Lagrangian  framework  with  a  maximum  CFL 
number  of  4.  The  total  error  in  the  mean-square  sense  ( Ej )  and  the  maximum  norm  (E,x)  is  plotted  in 
the  lower  left  of  each  panel,  the  maximum  and  minimum  values  of  the  numerical  solution  (which  can 
be  compared  with  the  correct  values  of  0  and  1)  are  plotted  in  the  lower  right.  Both  the  Eulerian  and 
semi-Lagrangian  methods  are  performing  better  than  the  WENO  method.  At  least  as  significant  are 
the  timing  comparisons  for  this  test,  which  show  that  the  WENO  method  is  a  factor  of  5  slower  than 
the  Eulerian  selectively  limited  PCM  method  and  a  factor  of  20  slower  than  the  semi-Lagrangian 
PCM.  This  method  is  slightly  faster  and  considerably  more  accurate  than  the  semi-Lagrangian  FCT 
method  being  incorporated  into  the  WRF  model(Skamarock,  2006).  The  four  leftmost  panels  in  Fig.  9 
all  show  results  from  simulations  at  relatively  coarse  resolution  (50x50  gridpoints).  At  finer  resolu¬ 
tions  of  100x100  (Fig.  9c)  or  200x200  (Fig.  9f),  the  selectively  limited  PCM  method  yields  considera¬ 
bly  more  accurate  results  and  converges  in  proportion  to  the  fifth  order  of  the  grid  spacing. 

IMPACT/APPLICATIONS 

COAMPS  is  the  Navy’s  operational  mesoscale  NWP  system  and  is  recognized  as  the  key  model  com¬ 
ponent  driving  a  variety  of  DoD  tactical  decision  aids.  Accurate  mesoscale  prediction  is  considered  an 
indispensable  capability  for  defense  and  civilian  applications.  Skillful  COAMPS  predictions  at  resolu¬ 
tions  less  than  1  km  will  establish  new  capabilities  for  the  support  of  the  warfighter  and  Sea  Power  2 1 . 
Operational  difficulties  with  weapon  systems  such  as  the  Joint  Standoff  Weapon  (JSOW)  have  been 
documented  in  regions  with  fine-scale  topography  due  to  low-level  wind  shear  and  turbulence.  Im¬ 
proved  high-resolution  predictive  capabilities  will  help  to  mitigate  these  problems  and  introduce  poten¬ 
tially  significant  cost  saving  measures  for  the  operational  application  of  JSOW.  The  capability  to  pre¬ 
dict  the  atmosphere  at  very  high  resolution  will  further  the  Navy  sea  strike  and  sea  shield  operations, 
provide  improved  representation  of  aerosol  transport,  and  will  lead  to  tactical  model  improvements. 
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This  research  will  establish  important  direction  for  the  development  of  the  Navy's  next  generation  mi¬ 
croscale  prediction  system.  Emergency  response  capabilities  and  Homeland  Security  issues  within  the 
DoD  and  elsewhere,  such  as  LLNL,  will  be  enhanced  with  the  new  modeling  capability. 

TRANSITIONS 

The  next  generation  COAMPS  system  will  transition  to  6.4  projects  within  PE  0603207N  (SPA WAR, 
PMW-180)  that  focus  on  the  transition  COAMPS  to  FNMOC.  The  improvements  to  the  COAMPS 
dynamical  core  have  been  transitioned  to  the  SPAWAR  6.4  project  and  subsequently  to  operations  as  a 
result  of  the  marked  improvement  in  the  geopotential  height  bias  statistics. 

RELATED  PROJECTS 

COAMPS  will  be  used  in  related  6. 1  projects  within  PE  0601 153N  that  include  studies  of  air-ocean 
coupling,  boundary  layer  studies,  and  topographic  flows,  and  in  related  6.2  projects  within  PE 
0602435N  that  focus  on  the  development  of  the  atmospheric  components  (QC,  analysis,  initialization, 
and  forecast  model)  of  COAMPS. 
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Cloud  Droplet  Number  Density  Cloud  Droplet  Number  Flux 


Figure  1.  The  model  simulation  of  the  cloud  droplet  density  (left  panel)  and  the  turbulent  flux 
(right  panel)  of  the  droplet  density  in  a  DYCOMS II  observation  case.  Blue  lines  represent 
the  result  from  the  traditional  “down-gradient”  parameterization;  black  lines  that  from 
the  new  parameterization.  The  circles  are  the  leg  averages  of  the  aircraft  data  taken  in 
DYCOMS  II  flight  1;  and  horizontal  bar  is  the  scatter  of  the  data. 


Cloud  Droplet  Number  Density  Flux 


Figure  2.  Comparison  of  the  flux  between  COAMPS-LES  results 
(left  panel)  and  the  new  parameterization  (right  panel). 
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V 


Figure  3.  Mean  sea  level  pressure  bias  as  a  function  of  forecast  hour  for  Feb.  1-14,  2007 for  the 
control  (red)  and  new  lower  boundary  condition  for  vertical  velocity  (blue)  forecasts.  Only  the  fore¬ 
casts  initialized  at  00  UTC  were  included  in  this  figure. 
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Figure  4.  Grid  elements  (thick  black  lines)  and  computational  grid  (thin  blue  lines). 


8 


-120  -90  -60  -30  0  30  60  90  120 

Distance  (km) 


Figure  5.  Perturbation  vertical  velocity  (m  s'1)  after  24  h  of  integration. 
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Figure  6.  Collisional  efficiencies  of  hexagonal  ice  plates  colliding  with  supercooled  droplets 
as  (a)  computed  numerically  by  Wang  and  Ji  (2000)  and  (b)  as  parameterized  for  use 
in  COAMPS.  The  parameterization  shows  a  plot  of  the  collision  efficiency  for  a  single  crystal 
with  a  radius  of  800  pm.  The  red  bold  line  shows  a  extension  of  the  Wang  and  Ji  results  for 
droplets  greater  then  the  critical  droplet  radius  (as  determined  by  the  point  in  the  droplet 
spectrum  where  the  droplet  fall  speed  equals  that  of  the  ice  crystal). 
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Rain  Collecting  Snow 


Rain  Mixing  Ratio  (g/kg) 


Figure  7.  A  plot  of  the  fractional  depletion  of  snow  during  a  time  step  for  the  given  rain  water  mix¬ 
ing  ratio  values  and  a  model  time  step  of  6  s.  The  snow  mixing  ratio  was  set  to  a  fixed  value  of  1 
g/kg.  The  symbols  A,B,C,D  refer  to  the  various  sensitivity  runs  as  described  in  the  text. 


Figure  8.  Illustration  of  the  basic  behavior  of  the  exact  solution. 
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Figure  9.  Numerical  results  generated  by:  (a)  the  unlimited  PCM  method,  (b)  the  selectively  limited 
PCM  method  with  positivity  preservation,  (d)  the  standard  WENO  method,  and  (e)  the  selectively 
limited  PCM  method  using  a  semi-Lagrangian  framework  with  a  maximum  CFL  number  of  4.  The 
total  error  in  the  mean-square  sense  (E2)  and  the  maximum  norm  (EJ  is  plotted  in  the  lower  left  of 
each  panel,  the  maximum  and  minimum  values  of  the  numerical  solution  (which  can  be  compared 
with  the  correct  values  of  0  and  1)  are  plotted  in  the  lower  right.  Finer  resolutions  of 100x100  (Fig. 

9c)  or  200x200  (Fig.  9f)  are  shown  as  well. 


11 


